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The algorithm for the DPD fluid, the dynamics of which is conceptually a combination of molecular dynamics, 

fj Brownian dynamics and lattice gas automata, is designed for simulating rheological properties of complex fluids on 

1^ hydrodynamic time scales. This paper calculates the equilibrium and transport properties (viscosity, self-diffusion) 

C of the thermostated DPD fluid explicitly in terms of the system parameters. It is demonstrated that temperature 

I ' gradients cannot exist, and that there is therefore no heat conductivity. 

jrt I Starting from the A^-particle Fokker- Planck, or Kramers' equation, we prove an 7J-theorem for the free energy, 

*^ obtain hydrodynamic equations, and derive a nonlinear kinetic equation (the Fokker-Planck-Boltzmann equation) 
for the single particle distribution function. This kinetic equation is solved by the Chapman-Enskog method. The 
C^ , analytic results are compared with numerical simulations. 
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I. INTRODUCTION 

The interest of the last decade in dynamical and rheological properties of complex fluids has seen the introduction 
of several new numerical methods for carrying out computer simulations on hydrodynamic time scales, the simulation 
of which using molecular dynamics often results in intensive computational demands. These new techniques include: 
(i) lattice gas cellular automata (LGCA) ||;gl; (ii) lattice Boltzmann equation (LBE) and (iii) dissipative particle 
dynamics (DPD). 

The last method was introduced by Hoogerbrugge and Koelman y], and was modified by Espahol and Warren ^ 
to ensure a proper thermal equilibrium state. The primary goal of this paper is a theoretical analysis and explicit 
calculation of transport and thermodynamic properties in terms of model parameters. This is highly relevant in view 
of the current interest in applications of DPD to systems such as flows past complex objects Q, concentrated colloidal 
suspensions 0,0 , dilute polymer solutions ||,P and phase separation jl^l . 

The DPD algorithm models a fluid of N interacting particles out of equilibrium and conserves mass and momentum. 
Position and velocity variables are continuous, as in Molecular Dynamics (MD), but time is updated in discrete steps 
St, as in LGCA and LBE. The algorithm is a mixture of molecular dynamics, Brownian and Stokesian dynamics and 
LGCA's, with a collision and a propagation step. In the collision step each particle interacts with all the particles 
inside an action sphere of radius Rq through conservative forces F^j, dissipative forces F^j^y , which are proportional to 
both the stepsize 6t and a friction constant 7, and random forces F^^j, which supply the energy lost by the damping. 
Here i,j e {l,2,...iV} label the particles. In numerical simulations, this is implemented by simultaneously updating 
the velocities from their precoUision value v; to their postcoUision value v* according to the instantaneous forces 
exerted by all particles inside the action sphere. In the subsequent propagation step of fixed length St all particles 
move freely to their new positions ri{t + St) = ri{t) + \'*St. 

Usual forms of the conservative force mean that the particles may be considered as completely interpenetrable. 
These softer interactions have the computational advantage y,g] of allowing particle motion on the order of a mean 
free path Iq during each time step of fixed length 6t. This represents a substantial advantage over event driven MD 
algorithms for hard sphere fluids, where the length 5t of the free propagation interval is on average much shorter, 
especially at fluid densities. 

By ignoring some of the microscopic details of the interactions, which are presumably irrelevant for fluid dynamics, 
DPD has the advantages of LGCAs, but avoids the disadvantages of lacking Galilean invariance and of introducing 
spurious conservation laws. In fact, the "point particles" should not be considered as molecules in a fluid, but rather 
as clusters of particles that interact dissipatively |jj|] . The introduction of noise and dissipation represents a coarse- 
grained mesoscopic level of description and hydrodynamic behaviour is expected at much smaller particle numbers 
than in conventional MD. If ig — l/7"-^o denotes the characteristic kinetic timescale in DPD, with n — N/V the 
number density, d the number of dimensions and 7 the friction constant, then to is considered to be large compared 
to any molecular time scale. 

In this coarse-grained description the dominant interactions are the dissipative and random forces, whereas the 
conservative forces can be interpreted as weak forces of relatively long range and may be taken into account as a 
Vlasov mean field term in the kinetic equations. In addition, they can have the spurious effect of tending to force the 
DPD particles into "colloidal crystal" configurations, unless friction and noise are sufficiently large to prevent cooling 
into a lattice configuration BpUll]. In the second half of the paper, where we derive a kinetic equation for the single 
particle distribution function, the conservative force will be neglected. This corresponds to the strong damping limit 
(7 large). The random forces act effectively as repulsive forces to prevent collapse of DPD particles. 

A substantial contribution towards the understanding of the DPD fiuid was given by Espafiol and Warren p| , who 
derived a Fokker-Planck equation for the iV-particle distribution function in the limit of continuous time {St —>■ 0). 
These authors also modified the original algorithm by imposing the detailed balance conditions, which guarantee the 
existence of the proper thermal Gibbs' equilibrium, described by exp[— i?/0o] where H is the Hamiltonian of the corre- 
sponding conservative system and 9q — ksTo is the global equilibrium temperature. These results are briefly reviewed 
in Section || to establish the notation. Concerning the macroscopic evolution equations, Espafiol formally established 
PI the linearized Navier-Stokes equations and derived Green-Kubo formulae for the DPD transport coefficients using 
a Mori-Zwanzig projection operator technique. However, to date no quantitative evaluation of these formulae for 
DPD seems to exist. Hence, little is known explicitly about the approach to equilibrium, the validity of standard 
hydrodynamics (system size dependence, effects of generalised hydrodynamics) or about transport coefficients. For 
the transport coefficients, Hoogerbrugge and Koelman [Q have estimated the kinematic viscosity v = rj/p, where rj 
is the shear viscosity and p = nm is the mass density, as z^ ~ ^uRq with uRq ^ 1. This result has recently been 
extended in |l3,n3[ to include the bulk viscosity by applying the "continuum approximation" to the discrete equations 



of motion for the DPD particles, following suggestions of Hoogerbrugge and Koelman M. In Section [II we show how 



the free energy of the DPD fiuid monotonically approaches its equilibrium value by proving an i?-theorem for the 



Fokker-Planck equation of Espanol and Warren, and we make the connection with the detailed balance conditions 
derived in |5[|. 

As a first step towards establishing the full nonlinear hydrodynaniic equations we derive in Section H^ the full 
macroscopic conservation laws for mass and momentum density, as well as the energy balance equation (details are 
given in Appendix A). The conceptual basis for the existence of hydrodynamic equations is the local equilibrium state, 
which in DPD is very different from that in a molecular fluid, because of the unusual role of the temperature. In 
Section ^ we study in a quantitative fashion the decay of the energy density e(r, t) and "kinetic" temperature 9{r, t) 
towards thermal equilibrium with global temperature ^Oi a-nd we assess in what sense and on what timescale the 
DPD fluid describes an isothermal fluid out of equilibrium. This is done on the basis of a nonlinear kinetic equation 
-referred to as Fokker-Planck-Boltzmann (FPB) equation- for the single particle distribution function /(x, i). It will 
be obtained from the first equation of the BBGKY hierarchy for the DPD fluid in combination with the molecular 
chaos assumption. 



By solving in Section VI the FPB equation in the hydrodynamic stage, using the Chapman-Enskog method, we 
derive the constitutive relations and the Navier-Stokes equation. This enables us to calculate in Section VII the 
transport coefficients of shear and bulk viscosity, as well as the self-diffusion coefficient. 

So far, we have not discussed the discrete time version of DPD, as implemented in actual simulations. They show 
a sensitive dependence of thermodynamic and transport properties on time step 5t [|l2|-[l4| . A promising step towards 
understanding the (5t-dependence was recently taken by Marsh and Yeomans |14|, who calculated the equilibrium 
temperature as a function of the step size, determined stability criteria for the step size, and validated their result 
by extensive numerical simulations. We shall not attempt to present here a systematic study of the different 0{5t)- 
corrections to equilibrium distributions and transport properties, but postpone this to a later publication. 



The paper ends in Section VIIl with comments on the most important results and future prospects for DPD. 



II. THE FOKKER-PLANCK FORMALISM 



The dynamics of a DPD system defines the time evolution of an A^-particle system, specified by a point F — {xj = 
(vi, Yi)\i = 1, 2 • • • N} in phase space, in terms of stochastic differential equations. For a theoretical description it is 
more convenient to consider the equivalent Fokker-Planck equation, derived by Espanol and Warren pi . 

To interpret the separate terms in the Fokker-Planck equation, it is instructive first to consider the analogous 
Kramers' equation for the probability P(v, r, t) of a single particle of mass to, having a phase description x = (v,r) 
at time t: 
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The three terms on the right can be interpreted as follows. The first term is an external conservative force F(r) = 
— Vy(r). The term involving the damping constant 7 corresponds to the Langevin force —7V and the diffusive term 
with diffusion coefficient ^a^ results from the random force a^ in the equivalent Langevin description, which reads: 



dr 

dt 

dv F 

-77 = 7V - 

at m 



< 



(2) 



where a£, is Gaussian white noise with amplitude a and < ^ >= and < ^{t)i{t') >— lS{t — t'), where I is a 
d-dimensional unit tensor. 

If we impose that the stationary solution of the Kramers' equation be the Gibbs' distribution: Peg ^ exp{— (imv^-l- 
V{r))/9o}, then the diffusion coefficient must satisfy the following Detailed Balance (DB) condition: 



a^ = 



27^0 



(3) 



where ^o = ksTo is the temperature in thermal equilibrium, measured in energy units. 

The full Fokker-Planck equation derived by Espafiol and Warren for the DPD system is a direct extension of the 
Kramers' equation to N interacting particles. The time evolution of the A^-particle distribution function P{T,t) is 
governed by: 



dtP = {Cc + Cj, + /:r) p, 



(4) 



where the Conservative, Dissipative and Random parts of the evolution operator are defined respectively as: 






V-- 9 1 V-- F(R„) / d d 



dvi 2 '^-^ m \ dvi 9v,- 
d 



^2 



'^D = Yl ^^^oiRrj) f Ry • ^ j (Ry • V, 



Cn = E T-«(^'^) (^^^ • 7^ ) ^^^ • ( 7^ - 7^ ) ■ ^'^ 



The summations run over all particles and the only difference to the original [^j is that the parameters 7 and a 
have been scaled by the mass m such that 7 has the dimensions of an inverse time. The three terms above are the 
TV-particle extensions of the three terms on the rhs of (|l|). 

• The conservative part Cc results from the additive and central interparticle interactions due to a potential, 
V = ^ X^i j=ii4>{Rij) where R^j — Yi — Tj is the relative position and a hat denotes a unit vector. It is the 
Liouville operator for the corresponding conservative system and in the limit of zero noise and friction, equation 
(g|) reduces to the Liouville equation. 

• The second term is analogous to the dissipative term in the Kramers' equation. It accounts for the Langevin 
damping force between the pair {ij), which is proportional to the friction constant 7 and to the component of 
the relvative velocity v^j along the line of centres R^j, and is of finite range. This last property is described by 
a positive weighting function WD{Rij) that is only non- vanishing inside an action sphere of finite radius _Ro- 

• The last term in (ph represents the random noise and should be compared with the diffusive term in (Q). 
The random force aS^ij between the pair (ij) is directed along Ry and is proportional to awR{Rij) where the 
weighting function wnlRij) is again only non- vanishing within a finite action sphere. 

The ranges of the conservative, dissipative and random forces may all be different as the model stands. Moreover, one 
of the essential properties of DPD is that its dynamics conserves total particle number, N, and the total momentum, 
P = J2i 'TiVi. Consequently < N >= J (ix/(x, t) and < P >= J dxmvf{x, t) are constants of the motion. The latter 
is always set equal to as the total system is assumed to be at rest. Here /(x, t) is the single particle distribution 
function. 

In addition, we want to emphasize that microscopic momentum conservation is an essential property of a fluid 
model if it is to have a momentum density p(r,t)u(r,i) that is slowly varying in space and time. In contrast, we 
note that the energy of the system is not strictly conserved under the DPD algorithm. The equations for the mass, 
momentum and energy density will be discussed more fully in Section fV[ In typical applications the conservative 
force may be set equal to zero. The parameters 7 and a satisfy the equation (H) and the weighting functions are 
chosen such that W]j{r) = w'^{r), which constitute in combination with (H) the detailed balance conditions for the 



DPD system, as will be discussed in Section |II|. The density is typically chosen such that there are 5 to 10 particles 
within an action sphere, which means that the instantaneous total force on any particle is small on average. 

Figure 1 shows an enlargement of part of configuration space, showing a sequence of 20 consecutive particle positions 
evolving from a randomly chosen initial configuration. The trajectories are relatively smooth, in contrast to the 
discontinuous paths in hard core interactions, illustrating that the resultant force on each particle is relatively small 
at this parameter setting. 
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FIG. 1. Typical evolution of a particle configuration over a period of 20 timesteps, showing that DPD interactions are 
"soft" as compared to hard core interparticle interactions. The circle with radius Ro — 4 indicates the range of interaction. 



III. AN H-THEOREM 

Consider the following functional of the A^-particle distribution function P{T,t): 



T[P]= fdTP{T,t){H + eolnP{T,t)} 



where ^o = rna'^ /2j and H is the Hamiltonian of the corresponding conservative system: 



^ = E r'^" + ^ = E J-^' + ^ E '^(^ 



2 ^ ^'"'"' 



(6) 



(7) 



V is the potential energy and (j){Rij ) is the pair interaction. The functional can be interpreted as a sort of free energy 
T = E — OqS, where E —< H > is the average total energy and S = — <\nP > yields the total entropy. The goal 
of this section is to show that .F is a Lyapunov functional with dtJ-^ < and to investigate the implications of this 
result for the equilibrium solution of the Fokker-Planck equation. 
The time derivative of (o) yields in combination with (0): 



dtT = f dr{H + 9a In P + Oo} (Cc + Co + U) P(r, t). 



(8) 



We observe that the third term inside {• • •} in M) vanishes due to total probability conservation. Then consider the 
contribution {dtT}c to (||) due to the Liouville operator Cc- Partial integration with respect to r,; and v^ yields 
directly 



{dtT}c = - fdT {PCcH + e^CcP} . 



(9) 



Here CcH ~ {H, H} = because the curly brackets represent Poisson brackets, as can easily be demonstrated. The 
second term in (0) reduces to surface terms in r,; and Vi and therefore vanishes too. 



Next, we combine the remaining terms in m) and perform a partial v^ integration, also symmetrizing the result 
over i and j. The final result is: 

X [wu{R^,)^., ■ Vy + wUR^/^^^^^, ' (^ " Q^)) ^- (10) 

Now we make the following observation. If we choose: 

woir) ^wji{r) =w{r), (11) 

where 'w{R) is an arbitrary positive function vanishing for r > Rq, then 9tJF < 0, as the rhs of dlG) can be cast into 
the form: 

dtT = -^m^JdrP J2 w{Rrj) {■ ■ •}' < 0, (12) 

where {• ■ •} is the same as in (p^. Note that the equality sign apphes if and only if P is the solution of ( p^ ) below. 
Consequently, the free energy-type function J^[P] is a monotonically decreasing function of time, until it reaches 
equilibrium where P = Peq which is simply the solution of {• • •}=0 for every pair (zj): 
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Changing variables to the relative velocities of the particles, it is easy to prove that the equilibrium distribution of 
the system is separable in the velocities, and has the general form: 

Pe9(r)=A(ri,---,rw)expJ-2i^^m(v, -Uo)4 , (14) 

where Uq is a constant independent of r^ and t. We will only consider macroscopic systems which are not in uniform 
motion at long times and consequently limit ourselves to Ug = 0. 

The function Peq{T) is also the stationary solution of the Fokker-Planck equation (El) if A(ri, • ■ • ,rjv) satisfies 
CqA = 0. This yields the Gibbs' distribution as the equilibrium solution: 

1 f H 



Pe,(r) = -exp|--|, (15) 

where H is the Hamiltonian (|7|) of the system and Z is a normalisation constant. We assume that Peq is uniquely 
determined by the requirements that it satisfies (ITsI) and (0). Consequently the DPD system will always reach the 
same equilibrium state if left undriven, independent of the volume and number of particles. The temperature of this 
equilibrium state has a value 0q = ma'^ /2j, which only depends on the parameters of the model. So DPD describes 
a system, thermostated at 9q and with a free energy J-[Peq] at equilibrium. Note that, in contrast to the iJ-theorem 
for the Boltzmann equation (see e.g. |2^), no molecular chaos approximation is required to derive this result for the 
DPD system. 

In their original discussion ||5|, Espahol and Warren imposed that the Gibbs' distribution be the stationary solution 
of the Fokker-Planck equation. The consequences of this requirement can be seen by inserting (|l^) into (Q) . It leads 
to the so-called detailed balance constraint 

2 

woir) = i^-^wl{r) = wl{r) = w{r). (16) 

Z7C/0 



Consequently, the constraint imposed in (|11| ) is the detailed balance constraint for DPD. The important result from 
the _ff-theorem is that it demonstrates that the Gibbs' distribution (ITq) is the inevitable equilibrium distribution. 
Throughout the rest of the paper we shall restrict ourselves to dealing exclusively with DPD systems that obey the 
DB condition (O). If the DB condition is violated, no H-theoTem can be derived and the Gibbs' distribution is not a 
stationary solution of the FP equation for DPD. In this case, the stationary state of the system does not correspond 



to thermal equilibrium but to some driven state, which will in general exhibit long range spatial correlations (see e.g. 

The original version of DPD, introduced by Hoogerbrugge and Koelman [Q violates the DB requirement (|l6|) and 
therefore its stationary distribution will not approach a Gibbs' state but may exhibit spatial correlations, e.g. algebraic 
correlations r^** where d is the number of dimensions, extending far beyond the ranges of the conservative, dissipative 
and stochastic forces. This absence of thermal equilibrium is likely to be the reason for difRculties and inconsistencies 
discussed in [Q-§[. 

IV. MACROSCOPIC CONSERVATION AND BALANCE EQUATIONS 

In the previous section we have established the existence of and approach to a thermal equilibrium state for the 
DPD system that obeys the detailed balance condition (|lq). In this section, we address the problem of how the 
quantities of macroscopic interest evolve in time towards the final equilibrium state, concentrating on the local mass 
density p(r, t) — mn(r, t), the local momentum density p(r, i)u(r, t) and the local energy density e(r, t). 

As discussed in section O the microscopic dynamics of DPD conserve mass and momentum, and the corresponding 
macroscopic densities obey local conservation laws. As the total energy is not conserved under DPD, the evolution 
equation for the macroscopic energy density does not have the form of a local conservation equation, but contains 
source and sink terms corresponding to the random and dissipative forces respectively. In the final equilibrium state, 
these will balance each other. 

Consider a general macroscopic quantity < A >, defined through 

(A) = J dTAiT)P{T,t). (17) 

Its time evolution can be obtained from the Fokker-Planck equation (0) combined with the detailed balance condition 
(|l6|), by multiplying the Fokker-Planck equation with A{T, t), integrating over all F-space and performing one or two 
partial integrations with respect to v^ and v^ . The result is the general rate of change equation: 

a. (A) = (e (- • I- + I • ^) ^) - 7 ( E -(%) {^. • V.} {r. • I-} a) 



E-(i?.){R.-^}{a.-(^-^)}A 



^ ( }_^ u, (i?., ) ^ R., . — M R.. • ( — - — ) M ) (18) 



for any dynamic variable A{T). 

Consider the conserved mass density p{r,t) — mn{r,t) and the momentum density p(r, t)u(r,t) defined through: 

n{v, t)^(y25ir-ri}\= J dv/(v, r, t) 

nu{r,t)^/J2^r,S{r-ri)\= /" rfv/(v, r, t)v. (19) 

It is convenient at this stage to introduce the single particle and pair distribution functions, defined as: 

/(x,t) = /(v,r,f) = /E'5(x-x,)\ 
/(2)(x,x',t) = / E S{^ - x,)<5(x' - x,)\ . (20) 

Application of (llS) to the conserved densities in dl3) yields the macroscopic conservation laws: 

dtp = -V • pu 
dtipu) = -V-{puu + U), (21) 



where V = d/dr and IT is the local pressure tensor or momentum flux density in the local rest frame of the fluid. The 
continuity equation has been derived by setting A — J2i H'"^ ~ ^i) i'^to ([l8|). The only non- vanishing term is the one 
containing {d/dri)S{r — ri) — — V(5(r — r.^), and the continuity equation follows at once. Derivation of the conservation 
equation for the momentum density proceeds along similar lines by choosing A — ^^77iv,;5(r — r^). Details of the 
latter derivation are given in the Appendix H, where it is shown that 

n(r, t) = nK(r, t) + nc(r, t) + nD(r, i), (22) 

with kinetic(K), coUisional transfer(C) and dissipativc(D) contributions: 

Hk = /"dvTOVV/(v,r,i) 

He = i /dvdv' /"dRRF(R)/'^(v,r,v',r',t) 

Uu = -^mj f dwdv' /"dRw(i?){R-(v-v')}RR7^^^(v,r,v',r',i), (23) 

—(2) 

with R = r — r'. The kinetic flux contains the so-called peculiar velocity, V = v — u(r,t), and / is the spatially 
averaged pair distribution function: 

/'\v,r,v',r,i)= / dA/(2)(v,r + AR,v',r + (A-l)R,i). (24) 

Jo 



The kinetic and coUisional transfer contributions to the momentum flux in (G2[) and ( |23D are present in any particle 
model with conservative forces. They are dominant in systems with sufficiently high density - dense gases and liquids 
- where the potential energy contributions are non-negligible with respect to the kinetic fluxes. 

The explicit form for these coUisional transfer contributions is given in the literature for several cases: smooth 
potentials ||I^, elastic hard spheres |20J, or inelastic hard spheres pl|. The dissipative contribution IId results from 
the Langevin-type damping forces between the particles. The random forces do not contribute to the momentum flux. 



The iZ-theorem, derived in section III, guarantees the approach to thermal equilibrium, where the distribution 
functions take the form: 

/ \ d/2 f 2 

I m \ mv 

/(x) ^no^o(v) = no(^^ ^"PTW 

/(2)(x,x')=n2^o(v)^o(v')5(|r-r'|), (25) 

where g{R) is the pair distribution function in thermal equilibrium, no = N/V is the number density, and <po(v) the 
Maxwellian velocity distribution. 

The sum of the kinetic and coUisional transfer contributions reduces to the equilibrium pressure and the dissipative 
contribution vanishes. Thus, n — pqI with po given by the virial theorem: 

Po = "0^0 - ^ I dRR^g{R), (26) 

where 

F(R)^-R^. (27) 

Away from global equilibrium, the pressure tensor II(r, t) will contain the local equilibrium pressure and terms 
involving the viscosities. However, before the Navier-Stokes equations, or more generally, the full set of hydrodynamic 
equations, can be derived the concept of local equilibrium - which forms the conceptual basis of slow hydrodynamic 
evolution - has to be re-examined, as the energy is no longer a conserved quantity. This can only be done after 
identifying the slow and fast relaxation modes in DPD, on the basis of a kinetic equation. This will be done in section 

Next we consider the energy density, defined as: 
e(r,t) = /^e,(v)<5(r-r,)\ 

= /"dvimv2/(x,t) + i f dvdV /"dR<?!>(i?)/(2)(v,r,v',r-R,i) , (28) 



where ei(v) = \raV^ + X^iati 't>{Rij) is the microscopic energy per particle. Use of the rate of change equation ( [T^ ) 
leads after some lengthy algebra to the energy balance equation, as discussed in Appendix]^ It reads: 

5te = -V-q + r. (29) 

Here the explicit form of the source term |g^ is: 

r(r, ^) = 7 ( E ^(^^^) {^0 - f {^'^ • (^' - ^^)}'} '^(r - ^0\ , (30) 

where the term proportional to ^o is a source resulting from the random force, and the term with the minus sign is 
a sink resulting from the Langevin-type damping forces. In global equilibrium, the source and sink terms balance 
one another and Fegm = 0. The heat current q, given explicitly in equations ( A12- A15) of Appendix H, contains the 



standard kinetic and coUisional transfer contributions due to conservative forces, as well as dissipative contributions 
analagous to IId in (|2^), and q vanishes in global equilibrium. 

If r would be set equal to zero, equation (E9h would have the generic form of the energy balance equation in 
ordinary hydrodynamics, where the heat current would contain a term proportional to the temperature gradient. As 
will become apparent in later sections, this is not the case in the DPD system. Although ( p9| ) looks like a macroscopic 
equation for the energy balance in the presence of sources and sinks, it looses its physical significance after a relaxation 
time to, in which e(r, t) — )■ WQn{Y,t)^ and (p9|) reduces to the continuity equation. This will be discussed in Section 



VI, below equation ( [4lD 

One may also derive a balance equation for the free energy density, which would be a local version of the iJ-theorem 
of Section III or of the corresponding one of Section ^ for the Fokker-Planck-Boltzmann equation. It would enable 
one to identify the irreversible entropy production. A similar balance equation for the entropy density in a dilute gas 
can be derived from the Boltzmann equation [ [l5[ . 

V. FOKKER-PLANCK-BOLTZMANN EQUATION 

In this section we derive an approximate kinetic equation, referred to as the Fokker-Planck-Boltzmann (FPB) 
equation, for the single particle distribution function /(x, t), which is based on the molecular chaos assumption and 
has a collision term which is quadratic in /(x, t). Moreover, from here on the conservative forces will be neglected, 
which corresponds to the strong damping limit (7 large). 

This section is organised as follows. We start by deriving the first equation of the BBGKY hierarchy, which relates 
dtf to the pair function /(^'(x,x', t). Then the molecular chaos assumption: 

/(2)(x,x',t)~/(x,t)/(x',t), (31) 

yields a closed equation, the FPB equation, which again satisfies an iJ-theorem. Next we analyse the local equilibrium 
solution of the kinetic equation, which provides the conceptual basis for the existence of hydrodynamic equations and 
transport coefficients, as well as the justification for solving this kinetic equation for finding the "normal solution" by 
means of the Chapman-Enskog method. 

The first equation of the BBGKY hierarchy can be derived directly by applying equation dl8) to the /i-space density: 

/(x)=E'^(^-^^)- (32) 

i 

Its average yields /(x, t) on account of ([20|). The resulting equation of motion is : 

«./^P,r,o{E|v. + i:o."(«.,){-v.4 + ^j^(^-^)}/:ft.,,a 



Mj-^Mj / , 



(33) 



where the (:) contraction of tensors is defined by A : B = 'Yliaa ^a/sB/sa, with a, /3 denoting Cartesian components of 
vectors or tensors. The equation can be further simplified to: 



dtf = -V • ( E '^i^i'^ -^i)) + Ijr • ( E "^(^ ^ yii)w{Rij)'R.i.j JRy • Vy i 



' 9v , 
7^0 52 



dwdw . 



E w{R,j)%j%j6{-y^ - X,) ) , (34) 



where V — d/dr. Performing the integrals over all variables except x^ and Xj leads to: 



dtf + yr-Vf = j dV 



dK RRu;(R) : <! t^ (v - v') 

' av 



m 9v9v 



/(^)(v,r,v',r-R,i). 



(35) 



This is the first equation of the BBGKY hierarchy with the Fokker-Planck equation (0) taking the place of the 
Liouville equation as the evolution equation. Under the molecular chaos approximation (p3l|) we have the following 
closed equation for the one particle distribution function: 



dtf + V. V/ = /(/) =^ dv' dR RRu>(i?)/(v', r - R, t) 



-- V- V +— — — - 

(7V 771 OVOV 



/(v,r,i). 



(36) 



The molecular chaos approximation is a mean-field approximation, which neglects dynamical correlations resulting 
from correlated multiple collisions taking place inside an action sphere. As we have set all conservative forces equal 
to zero, the molecular chaos assumption is exact in the global equilibrium state. Indeed, simulation results show that 
this is in fact an excellent approximation in the small-time step limit (as shown in Figure 2). 
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FIG. 2. Two particle distribution function, where the separation is measured in units _Ro- The system parameters in the 
simulation were taken: N = 2000 particles, friction constant 7 = 1, particle density n — 0.2, action sphere radius Ro = 4, and 
wuir) = 2{l-r/Ro). 



It can be shown in a similar fashion to section III that the functional, 



^ = y"dx|imv2+0oln/(x,i)|/(x,t). 



(37) 



satisfies an 7?-theorem {dtJ-^ < 0) where the equality only holds if /(x,i) is given by the equilibrium form no'^o(v) of 
(Eq), which establishes the existence of a unique global equilibrium state. 

The next problem is to solve the nonlinear FPB equation, and to analyse the approach to equilibrium using the 
Chapman- Enskog (CE) method. According to this method, one can distinguish two stages in the evolution of the 
single particle distribution function /(x, i): a rapid kinetic stage and a slow hydrodynamic stage |23|. 

In the kinetic stage, /(x, t) decays within a characteristic kinetic time ig to the so-called normal solution /(v|a(r, t)) 
which depends on space and time only through the first few moments a(r, t) — J (iva(v)/(x, t), the conserved densities, 



where a(v) — {1, v- • •} are the collisional invariants. In fluid systems the time to is the mean free time, whereas in a 
DPD system, to is estimated from (|3q ) as tg ~ l/(7ni?o). 

In the subsequent hydrodynamic stage, / depends only on space and time through its dependence on the conserved 
densities. In this stage, the solution /(v|a(r,t)) of the FPB equation can be determined perturbatively, / = /o + 
A*/i + ■ • •, as an expansion in powers of a small parameter, /i ~ lo^ i which measures the variation of the macroscopic 
parameters over a characteristic kinetic length scale, Iq ~ tov — {l/j)y^9o/m, where v ~ ^/9oJrn is a typical mean 
velocity. Therefore the /i-expansion is essentially an expansion in the small parameter I/7 (c.f. solution to Kramers' 
equation in |2J). 

In the remaining part of this section, we focus on determining the lowest order solution /o of (p6|), which is the 
local equilibrium distribution. We first observe that the Ihs of (pq) is of 0{fi), as Off is proportional to dfU ^ 0{ijl), 
and similarly for the gradient term. The rhs of ( J3^ ) is of 0(1). This requires that, to the dominant order in /i, / 
should satisfy /(/o) = + 0(/i). To determine the solution /o, we delocalise the collision operator /(/o) by replacing 
/o(v', r — R, i) on the rhs of ( |36| ) by /o(v', r, t) + 0(/i). If we denote the delocalised collision operator by Iq, then /o 
is the solution of: 

/o(/o) = 0. (38) 

Guided by the iJ-theorem, we assume the standard form for the local equilibrium distribution: 



d/2 



^"^^'")="(S) '"P 



m(v — u)^ 
29 



(39) 



where n, 9 and u are arbitrary functions of r and t. Substitution of (39) into (38) shows, however, that the above /o 
is only a solution if 

(40) 



where 9o is the constant model parameter introduced below equation (|6|), which equals the global equilibrium tem- 
perature. The parameters n{r,t) and u(r,t) in (BS) are chosen to be the fluid density and flow velocity. Hence, 

a(r,t) = /"dva(v)/(x,t) = /"dva(v)/o(x,i) , (41) 

where a(v) — {l,v} is a collisional invariant. 

This observation (|4^) has a profound consequence on the physical processes occurring in the DPD system, and 
makes it very different from standard fluids with energy conservation. In fluids, there is a fast kinetic relaxation to a 
local equilibrium state specified by n(r, t), 9{r, t) and u(r, t), and a subsequent slow hydrodynamic relaxation of these 
fields to global equilibrium. The DPD system distinguishes itself from standard fluids in the sense that there is a fast 
relaxation on a time scale to to a local equilibrium state (^)-(^), specified by n(r, t), u(r, t) and a spatially uniform 
and constant temperature 9q. The subsequent slow relaxation involves only the density n(r, t) and fiow velocity u(r, t). 

Consequently, a DPD system is not able to sustain a temperature gradient on hydrodynamic time scales; there is 
no heat current proportional to a temperature gradient; and there is no heat conductivity. Thus, the DPD system 
describes a thermostated or isothermal process at a fixed temperature 9o ■ It may only model physical systems where 
the temperature either relaxes very rapidly to an equilibrium value or where the temperature is irrelevant (an athcrmal 
process). 

It is worthwhile to explore the differences between DPD and standard fiuids somewhat further. Recalling that 
conservative forces have been set equal to zero in the present situation, permits us to write the energy density, 
e(r, i) = {d/2)n{r,t)9{r,t) in terms of a kinetic temperature 9{r,t). Clearly n{r,t) is a slowly changing variable, but 
what is the behaviour of 9{r, t) ? 

To answer this question, it is sufficient to consider only small deviations from global equilibrium, Sf = f — noVo(v), 
and to linearize the FPB equation around noipoi'v)- From that equation we shall derive how 69{r,t) — 9{r,t) — 6*0 
decays to zero. Linearisation of ( pq ) yields, after some algebra: 

dtSf + V • V<5/ = uo-^ ■ fv + ^^) 5f + ^ / dIlwiR)RR : u(r - R)v^o(v) , (42) 

av \ m ov J 9o J 



where we have used the relation 



nou(r,t) = I dw5f{-x,t)v , (43) 
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and introduced the coefficient, ujq = 1/to 

.o = ^/rfR^i?)-^M. (44) 

The equation with its non-local integral operator on the rhs of ( p^ would be the starting point for studying generalised 
hydrodynamics with wave-number dependent transport coefficients. Here, however, we shall only consider the decay 
of 

S0{v,t)^^jd^(^-eo)sf . (45) 

Then the rate of change of 5f can be calculated from (p2|), and yields: 

dtSe + V ■— f dv ( ^ - 0o] ^Sf ^ -~2ujQSe. (46) 

The second term on the Ihs of (Effl is typically an 0{fi) correction to the dominant decay terms. So (Eq) is a simple 
relaxation equation which shows explicitly that the kinetic temperature 0(jc,t) decays within the kinetic stage to the 
global temperature 6*0 with a relaxation time tg = {2u!o)~^ = ^Iq. 

The conclusion is that the energy density in the hydrodynamic stage, given by e{r,t) = {d/2)9on{r,t), is still a 
slow but not an independent variable. It is strictly proportional to the density. Moreover, we can conclude that the 
free-energy-type functions (o) and (pTI ) represent the actual free energy of the DPD system in the hydrodynamic stage. 

Similarly, we can determine the local equilibrium part of the pressure tensor (ESh in absence of conservative forces, 
by replacing / in IIk by its local equilibrium form /o and f^'^'> in IId by /o/o according to the stosszahlansatz [ pO[ . 
To zerotli order in n, the dissipative part IId vanishes, and the local equilibrium pressure is given by 

Ho = f dv mVV/o = n(r, t)9ol. (47) 

These results will be needed in the next section to solve the FPB equation to linear order in /x. 

VI. HYDRODYNAMIC STAGE 
A. Chapman-Enskog Method 



In Section IV we have derived the Fokker-Planck-Boltzmann equation for the DPD system and described the 



Chapman-Enskog method for obtaining its solution /(v|n, u) in the hydrodynamical stage. The method requires that 



the rhs and Ihs of the FPB-equation (36) are expanded in powers of /x ~ Iq^ , using the expansion 

/(v|n,u) = /o + M/i + ... (48) 

Every V is replaced by ^V and the derivative dtfo is eliminated using the macroscopic conservation laws. The lowest 
order solution, which is the local equilibrium distribution. 






-^(v-u(r,t))^ 



(49) 



has been determined in the previous section. 

To obtain /i we expand the FPB equation in powers of /z, yielding 

^ + Mv • V/o = /(/o) + Kdl/df)fji + ... (50) 

We start with the rhs of (pOh, which has been calculated exactly to C'(//)-terms included. One finds after some algebra 
that /(/o) = 0{fi^). In the previous section it has only been verified that I(/o) = 0{fi). The latter result is not 
sufficient, whereas the former is sufficient for our present purpose. In the remaining terms on the rhs we replace the 
collision operator / by its delocalised form /q, as defined below (pS]). The rhs of (|5^) then becomes: 

11 



W#)/o/i=-o^-(v+^^)A, (51) 



with Wo defined in ([l4[). 



To calculate the Ihs of (50) to 0{fi) we need the rate of change of n and u to lowest order in fi, which may be 



calculated from the conservation equations (|2l| ) with 11 replaced by its local equilibrium part IIo = ijOqI, calculated 



in (47), i.e 



dtn = -V • (nu) (52) 



dtu = — u ■ Vu Vn. 



They yield in combination with (Efl 



ith 



^+v-V/o = /o[J:D + JV-u] (53) 



JaM^) = ^ [va^Vp - ^6^pV' } (54) 



If 2 



and 



Note that the density and temperature gradients are absent on the rhs of (m) , in contrast to the traditional Chapman- 
Enskog result po| . 

For convenience of notation wc introduce the Fokker-Planck operator C and its adjoint £+, defined as 

/:=^-[V+^A] (56) 

£+ = r V + ^— 1.— 

to write the final equation for /i as: 

LOoCfi ^ fo[^ : D + JV ■ u]. (57) 

It is a second order PDE with an inhomogeneity on the rhs. We first construct a special solution by recalling that 
the Fokker-Planck operator £ can be mapped onto the Schrodinger equation for an isotropic d-dimensional harmonic 
oscillator [M. Its eigenfunctions are the tensor Hermite polynomials, usually called Sonine polynomials in a kinetic 
theory context, and the microscopic fluxes J and J are among them, i.e. 

r/oJ = -2/oJ ; r+J = -2J (58) 

CfoJ - -2foJ ; C+J - -2J. 



This can easily be verified. Combination of (57) and (|58|) yields the special solution 



/i--T^/o[J:D + JV-u]. (59) 

The general solution is obtained by adding an arbitrary linear combination of coUisional invariants a(V) — {1,V}, 
which are the solutions to the homogeneous equation £/oa(V) = 0. However, the constraint ( plj ) suppresses these 
terms and /i is the desired solution of the FPB equation to linear order in fi. 
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B. Navier-Stokes Equation 

The only slow macroscopic fields are the density n{r,t) and the flow velocity u(r, i), leading to the continuity 
equation and Navier-Stokes equation. The energy density in the hydrodynamic stage, e(r, t) = {d/2)n{r,t)9o, is not 
an independent variable. The energy balance equation derived in the Appendix |A| is only relevant in the kinetic stage, 
but has no physical significance in the hydrodynamical stage. 

The results for /o and /i are sufHcient to obtain the hydrodynamic equations to Navier-Stokes order and to obtain 
explicit expressions for the transport coefficients. The 0(/i)-correction /i contains only gradients of the flow field, 
Vu, but no gradients of the temperature. Therefore, there will be no heat current and vanishing heat conductivity. 
The C'(^)-terms in the pressure tensor 11 = XIq + /illi -I- ... will be proportional to Vu and we define the viscosities 
as the coefficients of proportionality through the constitutive relation 

Hi = -2t]D - CV ■ ul, (60) 

where 77 and ^ are respectively shear and bulk viscosity. 

Combining (60), (E^) and (pl|) then yields the Navier-Stokes equation for the DPD system 



dtipu) + V • (puu) = -OoVn + V • (27^0 + (V ■ ul). (61) 

The explicit expressions will be obtained in the next section. 

VII. TRANSPORT COEFFICIENTS 

A. Kinematic viscosities tjk and C^k 

There are two contributions to the pressure tensor: the kinetic part IIk and the dissipative part IId, defined in 
(23) and two corresponding viscosities. The kinetic part depends only on / = /o + /i/i, which are given in ( [49[ ) and 
(59). Then IIk becomes 

nK = n6'oI + MnK,i + -- (62) 

where /i is a formal expansion parameter that will be set equal to unity at the end of the calculations and 

Hk,! = J dvmVV/i =eafdw {J(V) + ^(V)I} A. (63) 

Here the dyadic mW has been split up into a traceless tensor, 6qJ, and a term 0oJ7I, proportional to the unit tensor 
and we have used the relation J dvfi — (see ([4l|)). Inserting the explicit solution ( p9|) into ( p^ ) allows us to write 
(p2) in the form 

UK,l = -^<^\J>■.D-^<J\J>V■nl, (64) 

ZujQ Zluq 



where we have used the relation /o = nfa{V) (see (25)) and introduced the inner product 



< A\B >^ / dvipo{V)A(V)B{V). (65) 

Moreover a crossproduct of a traceless tensor and a scalar vanishes , i.e. < J| ^7 >= 0. The product < J^\J^ > involves 
a simple Gaussian integral and yields 

<J\J>^^. (66) 

The fourth rank tensor < J|J > is isotropic, traceless and symmetric, which implies the general form 

< Ja^IJi^ > = ("p") fdv^V) (vc^Vp - ^S^pvA V-^Vs (67) 



C 



s.,s,..s.,,.-ls.,s,, 
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By taking double contractions and evaluating Gaussian integrals the constant C comes out to be equal to 1 and the 
(:) product in ( |64|) yields 

<J|J>:D = 2D. (68) 

Combination of (^^, (^6|) and ( pSf ) finally yields 

uOq n9o 
IIk 1 = D — V • ul. (69) 

ujq acjQ 



Comparison with the constitutive relation (|60|) enables us to identify the coefficients as the kinetic parts of the 
viscosities 

VK = T^— = TTTT ''K = -1— = -p-r , (70) 

where the definition (p^ ) of ujq has been used. We note that the kinetic part is inversely proportional to 7 and has 
been for the first time explicitly calculated. 

B. Dissipative viscosities t^d and C^o 

This section deals with the dissipative part 11 of the pressure tensor in ( p3| ) , which depends on the pair distribution 
function f^'^'. This function has a local equilibrium part /q and a part /i/| , linear in the gradients. We start with 
the first part. 

In order to make a direct comparison with the work of Espaiiol m we retain the conservative forces, for the time 
being. Then, the local equilibrium pair function has the form 

/(')(x,x') = /o(x)/o(x').9o(|x-x'|), (71) 

where goiR) is the spatial correlation function in local equilibrium. Only at the end of the calculation we will set the 
conservative forces equal to zero, so that go{R) = 1. 
Substitution of ( |7l| ) into (23) yields then 

Hd = -57m / dRgo(i?)w(i?)RR (R • [u(r) - u(r - R)]) ?i(r)n(r - R) (72) 

~ -^jmn^ / dRi?2.go(^)w(i?)RRRR : Vu , 

where [...] has been expanded to linear order in the gradients. Calculation of the completely symmetric isotropic 
fourth rank tensor proceeds as in ( p7|) with the result 

/ dRR^ ga(R)w{R)RaRf3R',R& = , ° [5ap5^s + 5a&5p^ + Sa^Sps] , (73) 

where 

[R^wgo]= I dRR^go{R)w{R). (74) 

We note that the definition of IId in (p3) contains / rather than /'■^•'(v, r, v', r — R, t). One easily verifies that the 
spatial averaging, denoted by the overline, makes no difference to linear order in the gradients. The final result for 
the dissipative part then becomes 

^° - did + 2) ^ 2d^ ^ ■ "^- ^^^> 

With the help of (|o|) the coefficients can be identified as the contributions to the viscosities due to the dissipative 
forces, i.e. 
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VB 



771771^ [i?^wgo] 



Cd 



771777^ [i?^w(7o] 



(76) 



2d{d + 2) ^" 2(P 

The local equilibrium contribution ( |75| ) to the dissipative pressure tensor turns out to be the dominant contribution 
to the viscosity of a DPD fluid, for large values of 717, as illustrated in Figure 3 and confirmed by numerical simulation 
in |4|, |l2|]. We also want to point out that an 0(/x)-contribution to the pressure tensor, calculated in local equilibrium 
as in (|72|), is not a novelty of this paper, but also occurs in all systems with impulsive (hard core) interactions that 
are not strictly local. For instance, consider the coUisional transfer contribution analogous to lie for elastic hard 
spheres, where F = — Vt/f in (E3h is ill-defined. This term is calculated in sections 16.4 and 16.5 of fpO|, where its local 



equilibrium contribution yields Tyns 



|Chs = t^ witn w 



with 



defined in Eq. (16.5.7) of ||20(]. These contributions 



,N 



n 
in real fiuids are the direct counterparts of 77D = |Cd "^ n'^ in DPD. 

We return to the DPD system without conservative forces, where the Gibbs' distribution ( [l5|) reduces to YiiLi fo{vi), 
and where the spatial correlations are absent, i.e. goiR) — 1- This equality is also required here for consistency with 
the molecular chaos approximation (pl|), used in section [V| and subsequent ones. 

So far, we have only considered local equilibrium contributions to IIo in (|7^). To obtain the complete contribution, 
consistent with the molecular chaos assumption, we substitute ( pl[ ) into (23) and use the definitions (|l9|). Surprisingly, 
the results (f2|) are recovered, showing that (|7^)-(|7l) give the full contribution of IId to the Navier-Stokes equation, 
at least within the molecular chaos assumption. 

To facilitate the comparison with the original predictions of p|p^p^ , we set go = 1 in (|7^) and introduce 



< i?2 >^= [R^w]/[wi 

where [a] denotes the spatial average introduced in (|44|), so that < R^ >^ 
The final result for the dissipative part of the viscosities is then 



R-o- 



m 



Cd 



777777^ < R^ >y 



\w\ 



2d{d + 2) 
777777^ < i?^ >„ 



2d2 



ujQtlneo/2[d + 2) 



= ujotl^n0o/2d, 



(77) 



(78) 



where wq = 1/io = "yn[w]/d is the characteristic relaxation rate introduced in (Q), and t^, defined through i^ = 
< R^ >w v'^, is the average trasversal time of an action sphere with v = {60/my^^ the thermal velocity. These results 
are in fact the theoretical predictions for the total shear and bulk viscosity of the DPD fluid, as obtained in [u2 13 
on the basis of the "continuum approximation" to the equations of motion of the DPD particles. In the present 
context of non-equilibrium statistical mechanics and kinetic theory, these contributions have been identified as the 
local equilibrium contributions to the transport coefficients in order to make the connection with Hoogerbrugge and 
Koelman's expression for the kinematic viscosity: 



v = v/p^ 



uj < R^ >w 
2d{d + 2)St 



(79) 



with their friction constant lu — ^St, proportional to St, as a proper friction should. Moreover, we recall that the 
range function w{R) in [Bj is normalised as. 



Tl\W\ 



dliwiR) = 1. 



(80) 



So, the results ([78) and (79) are identical. Hoogerbrugge and Koelman have also shown that the viscosity found in 
their numerical simulations approaches (^8|) and ( [79| ) for large 717. Simulations carried out with the modified DPD 
algorithm show the same properties |l2| . 

We conclude this subsection by listing the full results ( [70| ) and (^8|) for the shear and bulk viscosity in a DPD fiuid 
with continuous time (5t -^ 0): 



■n^m 



m = -iznOo 



C = Cd + Ck = ^ni 



d + 2 

^atl 



1 
1 



(81) 



and 



They involve the two intrinsic time scales of the DPD fluid: the characteristic kinetic time to = 1/ujq (see (jA 
the traversal time tyj of an action sphere, as defined below ( [78|), which is of order Rq/v. 

In the parameter range t^ > to the estimates 7]d and ^j of [n3,n3| dominate, and in the range t^ < to the kinematic 
viscosities do, as illustrated in Figure 3. 
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C. Numerical simulations 

As a simple test, the shear viscosity 77 of the DPD system was m.easured in two dimensions using a physical method. 
A linear velocity gradient was established between two moving plates and the force required to maintain this system 
was measured once equilibrium had been attained. 

By means of these simulations, we have measured the viscosity of the DPD fluid as a function of 717 at different 
temperatures ^o = ma'^ /2^. Results are shown in Figure 3 for a higher temperature to emphasize the importance of 
the kinematic contribution. At large 717 the measured viscosity approaches the theoretical prediction when the time 
step 5t is reduced. In this range of parameters, the viscosity is dominated by its dissipative part ^^ , corresponding to 
the original estimates of Hoogerbrugge and Koelman. At small 717 and high temperature Oo the viscosity is dominated 
by the kinetic viscosity. 

At small ri7 there are sizeable differences between predicted and simulated results, which do not decrease with 
decreasing time step size. The breakdown of the theory in this range of parameters has a fundamental reason. 
Inspection of the collision term on the rhs of (pq) or ( |5l| ) shows that with wg ~ n'y and loq9q ^ a^ small the typical 
size of the collision term ~ 1/to may not be large compared to the propagation terms on the rhs of (pq). Consequently, 
the Chapman-Enskog expansion will be poorly convergent or even divergent, because the kinetic and hydrodynamc 
time regimes are no longer well separated, or, equivalently, because the change of the macroscopic flow velocity over 
the characteristic kinetic length scale becomes large. To be consistent with the physical requirement of well separated 
time scales in this range of parameters, the imposed velocity gradients have to be reduced. 

Remaining differences between theory and simulations may be caused by a breakdown of the molecular chaos 
assumption, which neglects the dynamical correlations that may have to be taken into account. 
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FIG. 3. Kinematic viscosity v = r]/p against density n for dt — 0.05 and dt = 0.15. The system parameters in the 
simulations were taken : friction constant 7 = 1 and strength random force a = 1.5 for densities n = 0.025, 0.1, 0.3 and 0.4. 



D. Self-diffusion coefficient D 



The coefficient of self-diffusion can be obtained by considering a DPD fluid that is in equilibrium, except for the 
probability distribution fs of a tagged particle, labeled as i = 1. Following the arguments of Section M and choosing 
the /i-space density /s(x) — 5{x — xi), instead of (p^), one arrives at an equation similar to (p5|), with / and f^^'' 
replaced by fs and fs respectively, defined as 
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/,(x,t) = <<5(x-xi)> (82) 

/P(x,x',i) = < ^<5(x-xi)<5(x-x,) > . 

The molecular chaos assumption ( pi] ) now takes the form 

/p)(x,x',t) = n^<,(«')/s(x,t), (83) 

as the fluid particles are in thermal equilibrium with the Maxwellian (po{v) defined in (Eq), and the resulting FPB 
equation is linear, i.e. 

5t/.+v.V/.=c.o|--[v+^|-]/,. (84) 

It is identical to the Kramers' equation (|y) with F(r) = 0. 
The continuity equation takes the form 

5tc(r,t) + V-j(r,t) = 0, (85) 

with tagged particle density and current defined as 



c(r,i)=y dv/,(x,i), (86) 

j(v,i) = y"dvv/4x,t). 

Application of the Chapman-Enskog method to ( p4| ) yields the "local equilibrium" distribution function fso = 
c{r,t)ipo{v) and following equation for /si = fs — fso, 

B OB 

(Po{v)v ■ Vc = UJQTT- ■ (v + — -^)/sl = i^QjCfsl- (87) 

av m av 

As V0O at the Ihs is again an eigenfunction of C with eigenvalue —1, we find 

/si = - — 0o(«)v-Vc. (88) 

The coefficient of self- diffusion D, defined through the constitutive equation: 

j = / dvv/i = -DVc, (89) 

becomes for the DPD fluid: 

where p — mn is the mass density of the fluid. The above result is new. There is only a kinetic contribution and no 
dissipative one. 

VIII. CONCLUSIONS AND PROSPECTS 

The main results of this paper are the derivation and solution of the Fokker-Planck-Boltzmann (FPB) equation 
for the DPD fluid, providing explicit results for the thermodynamic and transport properties in terms of the system 
parameters: density n, friction constant 7, temperature Oq — ma'^ /2^ and range function 'w{R) with range i?o- There 
are two intrinsic time scales: the kinetic relaxation time t^ ^ I/ujRq determined by the collision term, and the 
traversal time i^, ~ Rq/v of an action sphere, where v = y/Oo/rn is the average velocity. We highlight the most 
important new results and future prospects in a number of comments. 
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1. The DPD fluid for continuous time (step size Ji — > 0), described by the iV-particle Fokker-Planck equation of 
Espahol and Warren, obeys an i/-theorem for the free energy T . The indispensable role of the detailed balance 
condition in establishing such a theorem is demonstrated. It guarantees a monotonic approach of T towards a 
unique thermal equilibrium, described by the Gibbs' distribution with a temperature 0o = ma'^ /2'y. 

2. The local conservation laws for mass density p = nm and momentum density are the essential prerequisites for 
the validity of the Navier-Stokes equations. The temperature, however, plays a very peculiar role. On the one 
hand the detailed balance condition guarantees the existence of a well-defined thermal equilibrium with a global 
equilibrium 9q, in which energy is conserved on average. On the other hand, the local equilibrium state depends 
only on n(r, t) and u(r, t), but not on a local equilibrium temperature 0(r, i), which relaxes in a time to (kinetic 
stage) towards its uniform equilibrium value 9q. In the subsequent hydrodynamic stage the DPD fluid is not 
able to sustain a temperature gradient, there is no heat conduction and all processes occur isothermally 

3. In the coarse-grained mesoscopic interpretation of DPD particles as "lumps of fluids" , the microscopic conser- 
vative forces between the DPD particles are small compared to the mesoscopic friction and random noise (large 
7 limit), and have been neglected in deriving the FPB equation. At sufficiently low temperature conservative 
forces can have the effect of forcing the DPD particles into crystalline configurations. 

4. The FPB equation is derived from the first equation of the BBGKY-hierarchy for the distribution functions, 
obtained from the iV-particle Fokker-Planck equation, playing the role of the Liouville equation. In addition, 
the molecular chaos assumption, /'^^^(x, x') — /(x)/(x') has been used. 

5. The Chapman-Enskog solution to the continuous time FPB equation yields two types of contributions to the 
viscosities: (i) Dissipative parts 77D and Cd, accounting for the coUisional transfer through the nonlocal dissipative 
interactions. They are determined by the local equilibrium distribution, (ii) Kinetic parts rjK and Ck, coming 
from the collision operator and determined by the Chapman-Enskog solution of the FPB equation. If t^ > to the 
dissipative viscosities are dominant; if t^ < to the kinetic viscosities are dominant. In ||l3,n3[ the total viscosity 
is estimated by rjn and Cd, which is correct for t^ ^ to. The simulated results for the kinematic viscosity are in 
reasonably good agreement with the predictions within the assessed theoretical regions of validity of the theory. 
We also calculated the coefficient of self-diffusion, which only has a kinetic part. 

6. It is also of interest to consider the Green-Kubo formulae for the viscosities in a DPD fluid as derived in [p|, 
where the linear 7-dependence of the viscosity in the limit of large dissipation 7 is questioned. To make the 
connection we observe that 77D and 77c of Jp[ should be identifled respectively with ?7d and tjk of the present 
paper. The time correlation functions of [^for ryD and r/c are formally proportional to 7^ and 1 respectively. 
Both time integrals appearing in the Green-Kubo formulae extend over the characteristic kinetic time to ^ I/7. 
Gonsequently, tjd '^ 7 and 77c '^ I/7 in the limit of large 7, in complete agreement with the detailed calculation 
of the present paper. 

7. The validity of the kinetic transport coefficients 77K and Ck and the convergence of the Chapman-Enskog ex- 
pansion require that spatial variations {fi ~ lo^) are small over a characteristic kinetic lengthscale lo ~ vto ~ 
v/wfRo. The convergence of the gradient expansion in ( |7^ ) and the validity of the dissipative viscosities ?7d 
and Cd require in addition that spatial variations are small over the diameter of an action sphere, Rq ~ vtw 
Both criteria pose bounds on the shear rates, imposed in the simulations, as well as on the validity of the 
Chapman-Enskog expansion. 

8. An interesting extension of the present theory would be towards generalised hydrodynamics. Such a region exists 
if tiu ^ to or Ro ^ lo- Then, the hydrodynamic modes with wave numbers k in the range (27r/i?o, 2tt/Io) have k- 
dependent dissipative viscosities ?7d and Cd • They may be calculated by studying the eigenmodes of the linearised 
FPB equation ([42[). A similar wavevector-range to generalised hydrodynamics occurs in dense hard sphere 
fluids, where Iq is small compared to the hard sphere diameter Ro- Such theories have been used successfully 
to describe neutron scattering experiments on liquid argon and liquid sodium |25| ]. Generalised hydrodynamics 
in DPD might therefore be of interest in explaining light and neutron experiments on concentrated colloidal 
suspensions. 

9. The equilibrium properties (see Figure 2 and Urn) and transport coefficients of DPD (see [ p^lJlql ) depend 
sensitively on the step size St. The Fokker-Planck equation (Q), the detailed balance condition (|16[), the FPB 
equation (pq), t he hydrodynamic equations (O) and corresponding transport coeffcients in subsections 1,2, and 
3 of Section^n] only hold for the continuous time model {6t — > 0) . The only analytic study, available on DPD at 
finite 6t J14| , calculates the equilibrium temperature 0{St), and derives criteria, imposed on 6t, for the stability 
of the equilibrium distribution /o(x). 
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The most important open problem on DPD is a systematic analysis of all (^^-corrections to equilibrium and transport 
properties, such as an explanation of Figure 2 and 3, suggesting that the current form of the modified DPD algorithm 
for finite step size St does not obey the detailed balance conditions, which implies that its stationary state is not the 
thermal equilibrium state described by the Gibbs' distribution. 
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APPENDIX A: DETAILS OF THE MACROSCOPIC FLOW EQUATIONS 

1. Momentum Conservation equation 



Inserting A = J2i m-Vi(5(r — rj in (18) yields directly: 



dtipu) = - V • / ^ mv,v,5(r " r,) \ + / ^ F,S{r - r,) \ - m7 / ^ w{R,j) (r^j.v,^) R,j,5(r - r,) \ . (Al) 

The first term on the rhs, which will be called rhsl, is transformed to the local rest frame of the fluid by introducing 
peculiar velocities, V^ = Vi — u(ri, t) and yields: 

rhsl == -V- (puu + Hk), (A2) 

where 

nK = femV,V,,5(r-r,)\ (A3) 

is the kinetic part of the pressure tensor, as listed in (|2^). The second term on the rhs of ( |Al[ ), which will be called 
rhs2, involves the conservative interparticle forces F,; = '^.j-a F(R-ij)- Symmetrising over i and j yields then: 



rhs2 - /^ E F(R..) [-^(r - r.) - <5(r - r,)]\ 

= -V • / i E F(R,,)Ry / dX S{r - r, + AR^ 



= -V • He. (A4) 

Here we have used the identity 

6{v ~ rO - ^(r - r,) = - J d\ —Sir - r, + ARy ) 

= -V • Ry / d\ 5{v - r, + ARy ). (A5) 



The third term on the rhs of (Al), referred to as rhs3, is due to dissipative particle interactions and can be treated 
in a sinnlar fashion. Symmetrising over i and j, and replacing ^(r — r^) by (1/2) [(5(r — r^) — S{r — Vj)] we obtain: 

rhs3 == ^ • ( y X! lw{Ri])'Rijtlij f Rij.Vy j / d\ S{y - r^ + AR^) ) 

\ i,i^i ° / 

= -V-nD. (A6) 
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The results ( k4 ) and ( [A^ ) are of the same general form, and can be expressed using the pair distribution function 
@ as: 

/^ A(R,j,v„v,)5(r-r, + ARy)\ = f dv f dV /'dRA(R, v, v')/''^(v,r + AR, v',r + (A - 1)R, t). (A7) 



where (A4), ([Aq ) and (A7) yield respectively He and IId as listed in (|23| ) with / defined in (P4[). Combination of 
(Al), (A2), (|A4[) and (A6) gives the macroscopic equation for the momentum density, 



?(2) 



(A8) 



9t (pu) + V • {puu + Hk + He + Hd} = 0, 
as listed in equations ( pl| ) and (|2^) in the body of the paper. 

2. Energy balance equation 

We start with the kinetic energy density ex by setting A ~ ^^(l/2)TOvf(5(r — r,) in (JTs]). This yields, after some 
algebra: 



r - r,; 



9teK = -V 7^ \m,v\b{r ~ r,)\ + ( J] v, • F(R,,)<5( 

\ i I \hj^i I 

- m7 / ^ w{Rij) (Ry • Vy) (Ry • V,) 5){v - r,) \ + 76I0 / ^ w{R^j)5{v - v,)\ . (A9) 

By setting A= ^ ^,- ■ /^ (j){Rij)S{r — r^) we find similarly for the potential energy density, 

dte^ = -^-Uj2 v.^l^.jO'Jlr - r0\ - i / ^ v,,.F(R,,)5(r - r,)\ . (AlO) 

We sum (A9) and ([AlCJ ) to obtain the rate of change of the total energy density: 

e = ex + 60 = / ^ ei(v)5(r - r.;) \ , 



where ei(v) is the microscopic energy per particle: 

e,(v) = -TOv,2 + -^0(A 






(All) 



(A12) 



We denote the n-th term on the rhs of (A9) and (AlC) by (an) and (bn) respectively and get the following results: 



(al) + (bl) - - V ■ / ^ v,e,(v)5(r - r,) \ = - V • qx 
(a2) + (b2) = \{Y. (v. + V,) • n^^i)6{v - rA 

= - V • / i ^ RyF(R„-) • (v, + vj) f dX S{r - r, + AR,j) 



ijT^i 



^V-qc. 



(A13) 



The expression for [(a2)+(b2)] has been symmetrized over i and j and ( |A5| ) has been used. The term (a4) represents 
the energy source Fr caused by the random forces. In (a3) we split v^ into (l/2)vij + (l/2)(vi + v^). The first term 
containing v^ gives the energy sink Fd resulting from the damping forces. The second term containing Vij is again 
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symmetrized and com.bined with (A5) to give the dissipative part of the energy current —V ■ qo- Combination of 
these terms then gives: 



(a4) + (a3) = 7 ( E ^i^'^M^ - ^0 ) - ^^"7 ( E ^(^»^) (^^^ ' ^^^)' ^(^ - ^0 ) 

+V ■ / - E R-y w(i?y) {fi^J ■ vy) Ry • (v, + Vj) / d\ 5{v - r, + ARy ) ) + 7 ( E w^(%)'5(r - r») \ 

= rR-rD-V-qD. (A14) 

where Fr, Fd and qd are defined by the three precedi ng te rms respe ctively. 



To obtain the full energy balance equation we sum ( All ) to (A14) to obtain 



9*6 = -V • [qK + qc + Qd] + Tr, - Fd 

= -V-q + F. (A15) 
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